# R version 4.0.3 (2020-10-10)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19045)
# 
# Matrix products: default
# 
# locale:
#   [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
# [5] LC_TIME=German_Germany.1252    
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] labelled_2.13.0  kableExtra_1.4.0 sandwich_3.0-0   lmtest_0.9-38    zoo_1.8-8        texreg_1.37.5    ggeffects_1.0.1 
# [8] haven_2.5.4      forcats_0.5.1    stringr_1.4.0    dplyr_1.0.3      purrr_0.3.4      readr_1.4.0      tidyr_1.1.2     
# [15] tibble_3.0.5     ggplot2_3.4.4    tidyverse_1.3.0 
# 
# loaded via a namespace (and not attached):
#   [1] fs_1.5.0          lubridate_1.7.9.2 insight_0.19.1    httr_1.4.2        rprojroot_2.0.2   tools_4.0.3       backports_1.2.1  
# [8] R6_2.5.0          sjlabelled_1.2.0  DBI_1.1.1         colorspace_2.0-2  withr_3.0.0       tidyselect_1.1.0  curl_4.3         
# [15] compiler_4.0.3    cli_3.1.1         rvest_0.3.6       xml2_1.3.2        labeling_0.4.2    scales_1.3.0      systemfonts_1.0.4
# [22] digest_0.6.27     foreign_0.8-80    rmarkdown_2.6     svglite_2.1.0     rio_0.5.16        pkgconfig_2.0.3   htmltools_0.5.1  
# [29] dbplyr_2.1.0      fastmap_1.1.0     rlang_1.1.2       readxl_1.3.1      rstudioapi_0.13   shiny_1.6.0       farver_2.0.3     
# [36] generics_0.1.0    jsonlite_1.7.2    zip_2.1.1         car_3.0-10        magrittr_2.0.1    Rcpp_1.0.12       munsell_0.5.0    
# [43] abind_1.4-5       lifecycle_1.0.4   stringi_1.5.3     yaml_2.2.1        carData_3.0-4     MASS_7.3-53       grid_4.0.3       
# [50] promises_1.1.1    crayon_1.3.4      miniUI_0.1.1.1    lattice_0.20-41   hms_1.0.0         knitr_1.48        pillar_1.4.7     
# [57] reprex_1.0.0      glue_1.4.2        evaluate_1.0.0    data.table_1.13.6 modelr_0.1.8      vctrs_0.6.5       httpuv_1.5.5     
# [64] cellranger_1.1.0  gtable_0.3.0      assertthat_0.2.1  xfun_0.47         ggExtra_0.10.0    openxlsx_4.2.3    mime_0.9         
# [71] xtable_1.8-4      broom_0.7.4       later_1.1.0.1     viridisLite_0.4.0 ellipsis_0.3.2    here_1.0.1 

# Libraries
library(tidyverse)
library(haven)
library(ggeffects)
library(texreg)
library(lmtest)
library(sandwich)
library(kableExtra)
library(labelled)


# Study to Load
dat <- readRDS("ExpStudy_00_data.RDS")

# Info about German 2021 Election ======


# RILE values for German parties (2021 Election)
party_data <- data.frame(
  party = c("Die Linke", "SPD", "Greens", "FDP", "CDU/CSU", "AfD"),
  rile = c(-36.167, -24.673, -21.038, 0.266, 3.495, 26.048)
)

# Plot the data
ggplot(party_data, aes(x = rile, y = 0)) +
  geom_point(size = 5, color = "black", alpha=0.3) +
  geom_hline(yintercept = 0, color = "black") +
  geom_text(aes(label = party), vjust = -1, size = 3) +
  scale_x_continuous(breaks = seq(-40, 30, by = 20)) +  # Adjust x-axis scale
  labs(
    title = "Left-Right Placement of German Parties (2021 Election)",
    x = "RILE Index (Left–Right)",
    y = ""
  ) +
  theme_minimal() +
  theme(axis.text.y = element_blank(), 
        panel.grid.major.y = element_blank(),  # Remove major horizontal grid lines
        panel.grid.minor.y = element_blank(),
        #panel.grid.minor = element_blank(),
        axis.ticks.y = element_blank())

ggsave("FigureSM4.pdf",width=9, height=4)


# Create the coalition government data frame
coalition_table <- data.frame(
  Election_Year = c(1961, 1965, 1969, 1972, 1976, 1980, 1983, 1987, 
                    1990, 1994, 1998, 2002, 2005, 2009, 2013, 2017, 2021),
  Coalition_Government = c("CDU/CSU-FDP", "CDU/CSU-FDP", "SPD-FDP", "SPD-FDP", "SPD-FDP", "SPD-FDP", 
                           "CDU/CSU-FDP", "CDU/CSU-FDP", "CDU/CSU-FDP", "CDU/CSU-FDP", "SPD-Greens", 
                           "SPD-Greens", "Grand Coalition", "CDU/CSU-FDP", "Grand Coalition", 
                           "Grand Coalition", "Traffic Light Coalition (Ampel)"),
  Parties_Involved = c("CDU/CSU, FDP", "CDU/CSU, FDP", "SPD, FDP", "SPD, FDP", "SPD, FDP", "SPD, FDP", 
                       "CDU/CSU, FDP", "CDU/CSU, FDP", "CDU/CSU, FDP", "CDU/CSU, FDP", "SPD, Greens", 
                       "SPD, Greens", "CDU/CSU, SPD", "CDU/CSU, FDP", "CDU/CSU, SPD", "CDU/CSU, SPD", 
                       "SPD, Greens, FDP")
)

# Print the table using knitr::kable for better formatting
latex_table <- kable(
  coalition_table, 
  format = "latex", 
  booktabs = TRUE, 
  caption = "Coalition Governments in Germany Since 1961",
  col.names = c("Election Year", "Coalition Government", "Parties Involved")
)

# Print LaTeX code
cat(latex_table, file = "TableSM19.tex")

# Quota Calculations ============

# Age
dat_age <- dat %>% 
  group_by(age) %>%
  summarise("n" = n()) %>%
  mutate("prop" = n/sum(n)) %>%
  na.omit() %>%
  mutate(var = paste("Age",to_character(age))) %>%
  select(var, n, prop)

# Region
dat_region <- dat %>% 
  group_by(state) %>%
  summarise("n" = n()) %>%
  mutate("prop" = n/sum(n)) %>%
  mutate(var = to_character(state)) %>%
  select(var, n, prop)

# Gender
dat_gender <- dat %>% 
  group_by(gender) %>%
  summarise("n" = n()) %>%
  mutate("prop" = n/sum(n))  %>%
  mutate(var = to_character(gender)) %>%
  select(var, n, prop)

# Bind Rows
output <- bind_rows(dat_gender,dat_region,dat_age) %>%
  na.omit() %>%
  kable(.,"latex",booktabs = T, linesep = "",digits = 2) %>%
  kable_styling(latex_options = "striped",stripe_index = c(1:3,4:(19),20:24)) %>%
  pack_rows("Gender", 1, 3) %>%
  pack_rows("Region", 4, 19) %>%
  pack_rows("Age", 20, 25)

cat(output, file = "TableSM20.tex")

# PAP Study I ===============

# Calculate mean and variance

## Computation of mean of CDU-specific coalition lottery (condition 2)
dat <- dat %>% 
  mutate(meanA_1 = as.numeric(Z_Y),
         meanA_2 = (as.numeric(Z_X)*as.numeric(p_2X)+as.numeric(Z_Y)*as.numeric(p_2Y)+as.numeric(Z_Z)*as.numeric(p_2Z))/100,
         meanA_3 = (as.numeric(Z_X)*as.numeric(p_3X)+as.numeric(Z_Y)*as.numeric(p_3Y)+as.numeric(Z_Z)*as.numeric(p_3Z))/100,
         varA_1 = 0,
         varA_2 = (as.numeric(p_2X)*(as.numeric(Z_X)-meanA_2)^2+ as.numeric(p_2Y)*(as.numeric(Z_Y)-meanA_2)^2 + as.numeric(p_2Z)*(as.numeric(Z_Z)-meanA_2)^2)/100,
         varA_3 = (as.numeric(p_3X)*(as.numeric(Z_X)-meanA_3)^2+ as.numeric(p_3Y)*(as.numeric(Z_Y)-meanA_3)^2 + as.numeric(p_3Z)*(as.numeric(Z_Z)-meanA_3)^2)/100,
         
         meanB_1 = as.numeric(Z_B),
         meanB_2 = (as.numeric(Z_A)*as.numeric(q_2A)+as.numeric(Z_B)*as.numeric(q_2B)+as.numeric(Z_C)*as.numeric(q_2C))/100,
         meanB_3 = (as.numeric(Z_A)*as.numeric(q_3A)+as.numeric(Z_B)*as.numeric(q_3B)+as.numeric(Z_C)*as.numeric(q_3C))/100,
         varB_1 = 0,
         varB_2 = (as.numeric(q_2A)*(as.numeric(Z_A)-meanB_2)^2+ as.numeric(q_2B)*(as.numeric(Z_B)-meanB_2)^2 + as.numeric(q_2C)*(as.numeric(Z_C)-meanB_2)^2)/100,
         varB_3 = (as.numeric(q_3A)*(as.numeric(Z_A)-meanB_3)^2+ as.numeric(q_3B)*(as.numeric(Z_B)-meanB_3)^2 + as.numeric(q_3C)*(as.numeric(Z_C)-meanB_3)^2)/100
) %>%
  mutate_at(vars(ends_with("Coal_text")), funs(gsub("B&uuml;ndnis 90/Gr&uuml;ne","Greens",gsub("<strong>|</strong>","",.)))) 
  
# Turn into long-format
dta_long <-  dat %>%
  mutate(id = 1:n()) %>%
  rename(#"vignetteA_0" = "ptv_union_1",
         "vignetteA_1" = "vignetteA1_1",
         "vignetteA_2" = "vignetteA2_1",
         "vignetteA_3" = "vignetteA3_1",
         #"vignetteB_0" = "ptv_greens_1",
         "vignetteB_1" = "vignetteB1_1",
         "vignetteB_2" = "vignetteB2_1",
         "vignetteB_3" = "vignetteB3_1",
         ) %>%
  select(starts_with("vignette"), id,starts_with("var"),starts_with("mean"), 
         ends_with("Coal_text"), starts_with("ptv"), risk_pref_1, id, coal_lik_union_1, 
         coal_lik_union_2, coal_lik_union_3, coal_lik_union_4,
         coal_lik_union_5, coal_lik_union_6, coal_lik_union_7, coal_lik_greens_1,
         coal_lik_greens_2, coal_lik_greens_3, coal_lik_greens_4, coal_lik_greens_5,
         p11,p12,p13,p14,p15,p16,p17,p21,p22,p23,p24,p25,p26,p27,p31,p32,p33,p34,p35,p36,p37,
         q11,q12,q13,q14,q15,q21,q22,q23,q24,q25,q31,q32,q33,q34,q35) %>%
  gather(var, ptv,starts_with("vignette"),starts_with("var"),starts_with("mean")) %>%
  separate(var, c("var","case")) %>%
  spread(var,ptv) %>%
  mutate(risk_pref = cut(risk_pref_1, 3),
         ptvA = cut(ptv_union_1,3),
         ptvB = cut(ptv_greens_1,3),
         sum_exp_union = exp(coal_lik_union_1) + exp(coal_lik_union_2) + exp(coal_lik_union_3) + exp(coal_lik_union_4) + exp(coal_lik_union_5) + exp(coal_lik_union_6) + exp(coal_lik_union_7),
         sum_exp_greens = exp(coal_lik_greens_1) + exp(coal_lik_greens_2) + exp(coal_lik_greens_3) + exp(coal_lik_greens_4) + exp(coal_lik_greens_5),
         coal_lik_union_1_norm = exp(coal_lik_union_1)/sum_exp_union,
         coal_lik_union_2_norm = exp(coal_lik_union_2)/sum_exp_union,
         coal_lik_union_3_norm = exp(coal_lik_union_3)/sum_exp_union,
         coal_lik_union_4_norm = exp(coal_lik_union_4)/sum_exp_union,
         coal_lik_union_5_norm = exp(coal_lik_union_5)/sum_exp_union,
         coal_lik_union_6_norm = exp(coal_lik_union_6)/sum_exp_union,
         coal_lik_union_7_norm = exp(coal_lik_union_7)/sum_exp_union,
         coal_lik_greens_1_norm = exp(coal_lik_greens_1)/sum_exp_greens,
         coal_lik_greens_2_norm = exp(coal_lik_greens_2)/sum_exp_greens,
         coal_lik_greens_3_norm = exp(coal_lik_greens_3)/sum_exp_greens,
         coal_lik_greens_4_norm = exp(coal_lik_greens_4)/sum_exp_greens,
         coal_lik_greens_5_norm = exp(coal_lik_greens_5)/sum_exp_greens,
         p11 = as.numeric(p11)/100,
         p12 = as.numeric(p12)/100,
         p13 = as.numeric(p13)/100,
         p14 = as.numeric(p14)/100,
         p15 = as.numeric(p15)/100,
         p16 = as.numeric(p16)/100,
         p17 = as.numeric(p17)/100,
         p21 = as.numeric(p21)/100,
         p22 = as.numeric(p22)/100,
         p23 = as.numeric(p23)/100,
         p24 = as.numeric(p24)/100,
         p25 = as.numeric(p25)/100,
         p26 = as.numeric(p26)/100,
         p27 = as.numeric(p27)/100,
         p31 = as.numeric(p31)/100,
         p32 = as.numeric(p32)/100,
         p33 = as.numeric(p33)/100,
         p34 = as.numeric(p34)/100,
         p35 = as.numeric(p35)/100,
         p36 = as.numeric(p36)/100,
         p37 = as.numeric(p37)/100,
         q11 = as.numeric(q11)/100,
         q12 = as.numeric(q12)/100,
         q13 = as.numeric(q13)/100,
         q14 = as.numeric(q14)/100,
         q15 = as.numeric(q15)/100,
         q21 = as.numeric(q21)/100,
         q22 = as.numeric(q22)/100,
         q23 = as.numeric(q23)/100,
         q24 = as.numeric(q24)/100,
         q25 = as.numeric(q25)/100,
         q31 = as.numeric(q31)/100,
         q32 = as.numeric(q32)/100,
         q33 = as.numeric(q33)/100,
         q34 = as.numeric(q34)/100,
         q35 = as.numeric(q35)/100,
         exp_diff_union_case0 = sqrt((coal_lik_union_1_norm-p11)^2+(coal_lik_union_2_norm-p12)^2+(coal_lik_union_3_norm-p13)^2+(coal_lik_union_4_norm-p14)^2+(coal_lik_union_5_norm-p15)^2+(coal_lik_union_6_norm-p16)^2+(coal_lik_union_7_norm-p17)^2),
         exp_diff_union_case1 = sqrt((coal_lik_union_1_norm-p21)^2+(coal_lik_union_2_norm-p22)^2+(coal_lik_union_3_norm-p23)^2+(coal_lik_union_4_norm-p24)^2+(coal_lik_union_5_norm-p25)^2+(coal_lik_union_6_norm-p26)^2+(coal_lik_union_7_norm-p27)^2),
         exp_diff_union_case2 = sqrt((coal_lik_union_1_norm-p31)^2+(coal_lik_union_2_norm-p32)^2+(coal_lik_union_3_norm-p33)^2+(coal_lik_union_4_norm-p34)^2+(coal_lik_union_5_norm-p35)^2+(coal_lik_union_6_norm-p36)^2+(coal_lik_union_7_norm-p37)^2),
         exp_diff_union = (exp_diff_union_case0+exp_diff_union_case1+exp_diff_union_case2)/3,
         exp_diff_union_cut = cut(exp_diff_union, 3),
         exp_diff_greens_case0 = sqrt((coal_lik_greens_1_norm-q11)^2+(coal_lik_greens_2_norm-q12)^2+(coal_lik_greens_3_norm-q13)^2+(coal_lik_greens_4_norm-q14)^2+(coal_lik_greens_5_norm-q15)^2),
         exp_diff_greens_case1 = sqrt((coal_lik_greens_1_norm-q21)^2+(coal_lik_greens_2_norm-q22)^2+(coal_lik_greens_3_norm-q23)^2+(coal_lik_greens_4_norm-q24)^2+(coal_lik_greens_5_norm-q25)^2),
         exp_diff_greens_case2 = sqrt((coal_lik_greens_1_norm-q31)^2+(coal_lik_greens_2_norm-q32)^2+(coal_lik_greens_3_norm-q33)^2+(coal_lik_greens_4_norm-q34)^2+(coal_lik_greens_5_norm-q35)^2),
         exp_diff_greens = (exp_diff_greens_case0+exp_diff_greens_case1+exp_diff_greens_case2)/3,
         exp_diff_greens_cut = cut(exp_diff_greens, 3)
         )

# Power Analysis ====

# Function to calculate MDE
calculate_mde <- function(n, alpha = 0.05, power = 0.9, sd = 1) {
  
  # Calculate the z-scores for the significance level and power
  z_alpha <- qnorm(1 - alpha / 2)
  z_beta <- qnorm(power)
  
  # Calculate the Minimal Detectable Effect
  mde <- (z_alpha + z_beta) * sd / sqrt(n)
  
  return(mde)
}


sigmaA <- sd(dta_long$vignetteA,na.rm=T)
calculate_mde(2*1577,sd=sigmaA) # Comparison prior and vignett 2 *n

sigmaB <- sd(dta_long$vignetteB,na.rm=T)
calculate_mde(2*1577,sd=sigmaB) # Comparison prior and vignett 2 *n


# Regression ====

# Model 1
  summary(m1A <- lm(vignetteA ~ case, as.data.frame(dta_long)))
  summary(m1B <- lm(vignetteB ~ case, as.data.frame(dta_long)))
  
  m1A_vcov <- vcovCL(m1A, type="HC1", cluster=dta_long$id)
  m1B_vcov <- vcovCL(m1B, type="HC1", cluster=dta_long$id)
  
  m1A_het <- coeftest(m1A,m1A_vcov)
  m1B_het <- coeftest(m1B,m1B_vcov)
  
  dfA <- ggpredict(m1A, terms = "case", 
                   vcov.fun = "vcovCL", 
                   vcov.type = "HC1",
                   vcov.args = list(cluster = dta_long$id)) %>% as.data.frame() %>% mutate(party = "CDU/CSU")
  dfB <- ggpredict(m1B, terms = "case", 
                   vcov.fun = "vcovCL", 
                   vcov.type = "HC1",
                   vcov.args = list(cluster = dta_long$id)) %>% as.data.frame() %>% mutate(party = "Greens")
  df <- bind_rows(dfA,dfB)
  
  df$group <- c("Certain", "Uncertain",
                "Very Uncertain")
  
  df_stars <- data.frame(
    est = c(coef(m1A),coef(m1B)),
    se = c(sqrt(diag(m1A_vcov)),
           sqrt(diag(m1B_vcov))),
    group = c("Intercept","Uncertain",
            "Very Uncertain"),
    party = c(rep("CDU/CSU",3),rep("Greens",3))
  )  %>%
    filter(group != "Intercept") %>%
    mutate(pval = pnorm(abs(est/se),lower.tail = F)*2 )%>%
    mutate(stars = case_when(
      pval < 0.001 ~ "***",
      pval < 0.01 ~ "**",
      pval < 0.05 ~ "*",
      TRUE ~ ""
    ))

  df <- left_join(df, df_stars, by=c("group","party" )) %>%
    as.data.frame()
  
  ggplot(df,aes(x=group,y=predicted,ymin=conf.low,ymax=conf.high)) + 
    geom_errorbar(width=.05) + theme_minimal() +
    geom_point() +
    geom_bar(stat = "identity", alpha=0.3) +
    facet_grid(~ party) +
    theme(text = element_text(size=16,family = "serif")) +
    ylab("Propensity to Vote") + xlab("Coalition Government Lottery") + 
    geom_text(aes(label = paste(ifelse(is.na(est),"",round(est, 2)),
                                ifelse(is.na(stars),"",stars)),y=conf.high),nudge_y = .2) 
  ggsave("Figure4.pdf", width = 8, height=4)
  

  # Fixed Effects
  summary(m1AFE <- lm(vignetteA ~ case + as.factor(id), as.data.frame(dta_long)))
  summary(m1BFE <- lm(vignetteB ~ case + as.factor(id), as.data.frame(dta_long)))

  m1AFE_vcov <- vcovCL(m1AFE, type="HC1", cluster=dta_long$id)
  m1BFE_vcov <- vcovCL(m1BFE, type="HC1", cluster=dta_long$id)
  
  m1AFE_he <- coeftest(m1AFE,m1AFE_vcov)
  m1BFE_he <- coeftest(m1BFE,m1BFE_vcov)
  
  # Scenarios
  summary(m1ASC <- lm(vignetteA ~ case + X_Coal_Text + Y_Coal_Text + Z_Coal_Text, as.data.frame(dta_long)))
  summary(m1BSC <- lm(vignetteB ~ case + A_Coal_Text + B_Coal_Text + C_Coal_Text, as.data.frame(dta_long)))
  
  m1ASC_vcov <- vcovCL(m1ASC, type="HC1", cluster=dta_long$id)
  m1BSC_vcov <- vcovCL(m1BSC, type="HC1", cluster=dta_long$id)
  
  m1ASC_he <- coeftest(m1ASC,m1ASC_vcov)
  m1BSC_he <- coeftest(m1BSC,m1BSC_vcov)
  
  # Model 2
  summary(m2A <- lm((vignetteA) ~ case*risk_pref, as.data.frame(dta_long)))
  summary(m2B <- lm((vignetteB) ~ case*risk_pref, as.data.frame(dta_long)))
  
  m2A_vcov <- vcovCL(m2A, type="HC1", cluster=dta_long$id)
  m2B_vcov <- vcovCL(m2B, type="HC1", cluster=dta_long$id)
  
  m2A_he <- coeftest(m2A, m2A_vcov)
  m2B_he <-  coeftest(m2B, m2B_vcov)
  
  dfA <- ggpredict(m2A, terms = c("case","risk_pref")) %>% as.data.frame() %>% mutate(party = "CDU/CSU")
  dfB <- ggpredict(m2B, terms = c("case","risk_pref")) %>% as.data.frame() %>% mutate(party = "Greens")
  df <- bind_rows(dfA,dfB)
  
  
  
  levels(df$x) <- c("Certain", "Uncertain","V. Uncertain")
  
  df$group <- factor(c("Low Risk Pref.", "Middle Risk Pref.",
                       "High Risk Pref."), levels=c("Low Risk Pref.", "Middle Risk Pref.",
                                                "High Risk Pref."))
  
  
  # Calculate Treatment Effects
  ateA_case2_ptvlow <- coef(m2A)["case2"]
  ateA_case2_ptvmid <- coef(m2A)["case2"] + coef(m2A)["case2:risk_pref(3.33,6.67]"]
  ateA_case2_ptvhgh <- coef(m2A)["case2"] + coef(m2A)["case2:risk_pref(6.67,10]"]
  
  ateA_case3_ptvlow <- coef(m2A)["case3"]
  ateA_case3_ptvmid <- coef(m2A)["case3"] + coef(m2A)["case3:risk_pref(3.33,6.67]"]
  ateA_case3_ptvhgh <- coef(m2A)["case3"] + coef(m2A)["case3:risk_pref(6.67,10]"]
  
  seA_case2_ptvlow <- sqrt(m2A_vcov["case2","case2"])
  seA_case2_ptvmid <- sqrt(m2A_vcov["case2","case2"] + 
                             m2A_vcov["case2:risk_pref(3.33,6.67]","case2:risk_pref(3.33,6.67]"] +
                             2 * m2A_vcov["case2","case2:risk_pref(3.33,6.67]"])
  
  seA_case2_ptvhgh <- sqrt(m2A_vcov["case2","case2"] + 
                             m2A_vcov["case2:risk_pref(6.67,10]","case2:risk_pref(6.67,10]"] +
                             2 * m2A_vcov["case2","case2:risk_pref(6.67,10]"])
  
  seA_case3_ptvlow <- sqrt(m2A_vcov["case3","case3"])
  seA_case3_ptvmid <- sqrt(m2A_vcov["case3","case3"] + 
                             m2A_vcov["case3:risk_pref(3.33,6.67]","case3:risk_pref(3.33,6.67]"] +
                             2 * m2A_vcov["case3","case3:risk_pref(3.33,6.67]"])
  seA_case3_ptvhgh <- sqrt(m2A_vcov["case3","case3"] + 
                             m2A_vcov["case3:risk_pref(6.67,10]","case3:risk_pref(6.67,10]"] +
                             2 * m2A_vcov["case3","case3:risk_pref(6.67,10]"])
  
  # Calcualte Treatment Effects
  ateB_case2_ptvlow <- coef(m2B)["case2"]
  ateB_case2_ptvmid <- coef(m2B)["case2"] + coef(m2B)["case2:risk_pref(3.33,6.67]"]
  ateB_case2_ptvhgh <- coef(m2B)["case2"] + coef(m2B)["case2:risk_pref(6.67,10]"]
  
  ateB_case3_ptvlow <- coef(m2B)["case3"]
  ateB_case3_ptvmid <- coef(m2B)["case3"] + coef(m2B)["case3:risk_pref(3.33,6.67]"]
  ateB_case3_ptvhgh <- coef(m2B)["case3"] + coef(m2B)["case3:risk_pref(6.67,10]"]
  
  seB_case2_ptvlow <- sqrt(m2B_vcov["case2","case2"])
  seB_case2_ptvmid <- sqrt(m2B_vcov["case2","case2"] + 
                             m2B_vcov["case2:risk_pref(3.33,6.67]","case2:risk_pref(3.33,6.67]"] +
                             2 * m2B_vcov["case2","case2:risk_pref(3.33,6.67]"])
  
  seB_case2_ptvhgh <- sqrt(m2B_vcov["case2","case2"] + 
                             m2B_vcov["case2:risk_pref(6.67,10]","case2:risk_pref(6.67,10]"] +
                             2 * m2B_vcov["case2","case2:risk_pref(6.67,10]"])
  
  seB_case3_ptvlow <- sqrt(m2B_vcov["case3","case3"])
  seB_case3_ptvmid <- sqrt(m2B_vcov["case3","case3"] + 
                             m2B_vcov["case3:risk_pref(3.33,6.67]","case3:risk_pref(3.33,6.67]"] +
                             2 * m2B_vcov["case3","case3:risk_pref(3.33,6.67]"])
  seB_case3_ptvhgh <- sqrt(m2B_vcov["case3","case3"] + 
                             m2B_vcov["case3:risk_pref(6.67,10]","case3:risk_pref(6.67,10]"] +
                             2 * m2B_vcov["case3","case3:risk_pref(6.67,10]"])
  df_stars <- data.frame(
    est = c(ateA_case2_ptvlow,ateA_case2_ptvmid,ateA_case2_ptvhgh,
            ateA_case3_ptvlow,ateA_case3_ptvmid,ateA_case3_ptvhgh,
            ateB_case2_ptvlow,ateB_case2_ptvmid,ateB_case2_ptvhgh,
            ateB_case3_ptvlow,ateB_case3_ptvmid,ateB_case3_ptvhgh),
    se = c(seA_case2_ptvlow,seA_case2_ptvmid,seA_case2_ptvhgh,
           seA_case3_ptvlow,seA_case3_ptvmid,seA_case3_ptvhgh,
           seB_case2_ptvlow,seB_case2_ptvmid,seB_case2_ptvhgh,
           seB_case3_ptvlow,seB_case3_ptvmid,seB_case3_ptvhgh),
    group = factor(c("Low Risk Pref.", "Middle Risk Pref.",
                     "High Risk Pref."), levels=c("Low Risk Pref.", "Middle Risk Pref.",
                                                  "High Risk Pref.")),
    party = c(rep("CDU/CSU",6),rep("Greens",6)),
    x=c(rep("Uncertain",3),rep("V. Uncertain",3))
  )  %>%
    mutate(pval = pnorm(abs(est/se),lower.tail = F)*2 )%>%
    mutate(stars = case_when(
      pval < 0.001 ~ "***",
      pval < 0.01 ~ "**",
      pval < 0.05 ~ "*",
      TRUE ~ ""
    )) 
  
  df <- left_join(df, df_stars, by=c("group","party","x" )) %>%
    as.data.frame()
  
  ggplot(df, aes(x=x,y=predicted,ymin=conf.low,ymax=conf.high)) + 
    geom_errorbar(width=.05) + theme_minimal() +
    geom_point() + 
    facet_grid(group ~ party) +
    geom_bar(stat = "identity", alpha=0.3) +
    ylab("Propensity to Vote") + xlab("Coalition Government Lottery") +
    theme(text = element_text(size=18,family = "serif")) + 
    geom_text(aes(label = paste(ifelse(is.na(est),"",round(est, 2)),
                                ifelse(is.na(stars),"",stars)),y=conf.high),nudge_y = .6) 
  ggsave("FigureSM5.pdf", width = 8, height=6)
  
  

  
# Model 3
  summary(m3A <- lm((vignetteA) ~ case*ptvA, as.data.frame(dta_long)))
  summary(m3B <- lm((vignetteB) ~ case*ptvB, as.data.frame(dta_long)))
  
  m3A_vcov <- vcovCL(m3A, type="HC1", cluster=dta_long$id)
  m3B_vcov <- vcovCL(m3B, type="HC1", cluster=dta_long$id)
  
  
  m3A_he <- coeftest(m3A, m3A_vcov)
  m3B_he <- coeftest(m3B, m3B_vcov)
  
  dfA <- ggpredict(m3A, terms = c("case","ptvA")) %>% as.data.frame() %>% mutate(party = "CDU/CSU")
  dfB <- ggpredict(m3B, terms = c("case","ptvB")) %>% as.data.frame() %>% mutate(party = "Greens")
  df <- bind_rows(dfA,dfB)
  
    levels(df$x) <- c("Certain", "Uncertain",
                    "V. Uncertain")


  df$group <- factor(c("Low PTV", "Middle PTV",
                       "High PTV"), levels=c("Low PTV", "Middle PTV",
                                             "High PTV"))
  
  # Calcualte Treatment Effects
  ateA_case2_ptvlow <- coef(m3A)["case2"]
  ateA_case2_ptvmid <- coef(m3A)["case2"] + coef(m3A)["case2:ptvA(3.33,6.67]"]
  ateA_case2_ptvhgh <- coef(m3A)["case2"] + coef(m3A)["case2:ptvA(6.67,10]"]
  
  ateA_case3_ptvlow <- coef(m3A)["case3"]
  ateA_case3_ptvmid <- coef(m3A)["case3"] + coef(m3A)["case3:ptvA(3.33,6.67]"]
  ateA_case3_ptvhgh <- coef(m3A)["case3"] + coef(m3A)["case3:ptvA(6.67,10]"]
  
  seA_case2_ptvlow <- sqrt(m3A_vcov["case2","case2"])
  seA_case2_ptvmid <- sqrt(m3A_vcov["case2","case2"] + 
                     m3A_vcov["case2:ptvA(3.33,6.67]","case2:ptvA(3.33,6.67]"] +
                     2 * m3A_vcov["case2","case2:ptvA(3.33,6.67]"])
  
  seA_case2_ptvhgh <- sqrt(m3A_vcov["case2","case2"] + 
                            m3A_vcov["case2:ptvA(6.67,10]","case2:ptvA(6.67,10]"] +
                            2 * m3A_vcov["case2","case2:ptvA(6.67,10]"])
  
  seA_case3_ptvlow <- sqrt(m3A_vcov["case3","case3"])
  seA_case3_ptvmid <- sqrt(m3A_vcov["case3","case3"] + 
                            m3A_vcov["case3:ptvA(3.33,6.67]","case3:ptvA(3.33,6.67]"] +
                            2 * m3A_vcov["case3","case3:ptvA(3.33,6.67]"])
  seA_case3_ptvhgh <- sqrt(m3A_vcov["case3","case3"] + 
                            m3A_vcov["case3:ptvA(6.67,10]","case3:ptvA(6.67,10]"] +
                            2 * m3A_vcov["case3","case3:ptvA(6.67,10]"])
  
  # Calcualte Treatment Effects
  ateB_case2_ptvlow <- coef(m3B)["case2"]
  ateB_case2_ptvmid <- coef(m3B)["case2"] + coef(m3B)["case2:ptvB(3.33,6.67]"]
  ateB_case2_ptvhgh <- coef(m3B)["case2"] + coef(m3B)["case2:ptvB(6.67,10]"]
  
  ateB_case3_ptvlow <- coef(m3B)["case3"]
  ateB_case3_ptvmid <- coef(m3B)["case3"] + coef(m3B)["case3:ptvB(3.33,6.67]"]
  ateB_case3_ptvhgh <- coef(m3B)["case3"] + coef(m3B)["case3:ptvB(6.67,10]"]
  
  seB_case2_ptvlow <- sqrt(m3B_vcov["case2","case2"])
  seB_case2_ptvmid <- sqrt(m3B_vcov["case2","case2"] + 
                             m3B_vcov["case2:ptvB(3.33,6.67]","case2:ptvB(3.33,6.67]"] +
                             2 * m3B_vcov["case2","case2:ptvB(3.33,6.67]"])
  
  seB_case2_ptvhgh <- sqrt(m3B_vcov["case2","case2"] + 
                             m3B_vcov["case2:ptvB(6.67,10]","case2:ptvB(6.67,10]"] +
                             2 * m3B_vcov["case2","case2:ptvB(6.67,10]"])
  
  seB_case3_ptvlow <- sqrt(m3B_vcov["case3","case3"])
  seB_case3_ptvmid <- sqrt(m3B_vcov["case3","case3"] + 
                             m3B_vcov["case3:ptvB(3.33,6.67]","case3:ptvB(3.33,6.67]"] +
                             2 * m3B_vcov["case3","case3:ptvB(3.33,6.67]"])
  seB_case3_ptvhgh <- sqrt(m3B_vcov["case3","case3"] + 
                             m3B_vcov["case3:ptvB(6.67,10]","case3:ptvB(6.67,10]"] +
                             2 * m3B_vcov["case3","case3:ptvB(6.67,10]"])
  df_stars <- data.frame(
    est = c(ateA_case2_ptvlow,ateA_case2_ptvmid,ateA_case2_ptvhgh,
            ateA_case3_ptvlow,ateA_case3_ptvmid,ateA_case3_ptvhgh,
            ateB_case2_ptvlow,ateB_case2_ptvmid,ateB_case2_ptvhgh,
            ateB_case3_ptvlow,ateB_case3_ptvmid,ateB_case3_ptvhgh),
    se = c(seA_case2_ptvlow,seA_case2_ptvmid,seA_case2_ptvhgh,
           seA_case3_ptvlow,seA_case3_ptvmid,seA_case3_ptvhgh,
           seB_case2_ptvlow,seB_case2_ptvmid,seB_case2_ptvhgh,
           seB_case3_ptvlow,seB_case3_ptvmid,seB_case3_ptvhgh),
    group = factor(c("Low PTV","Middle PTV","High PTV")),
    party = c(rep("CDU/CSU",6),rep("Greens",6)),
    x=c(rep("Uncertain",3),rep("V. Uncertain",3))
  )  %>%
    mutate(pval = pnorm(abs(est/se),lower.tail = F)*2 )%>%
    mutate(stars = case_when(
      pval < 0.001 ~ "***",
      pval < 0.01 ~ "**",
      pval < 0.05 ~ "*",
      TRUE ~ ""
    )) 
  
  df <- left_join(df, df_stars, by=c("group","party","x" )) %>%
    as.data.frame()
  
  
  ggplot(filter(df), 
         aes(x=x,y=predicted,ymin=conf.low,ymax=conf.high)) + 
    geom_errorbar(width=.05) + theme_minimal() +
    geom_point() +
    facet_grid(group ~ party) +
    geom_bar(stat = "identity", alpha=0.3) +
    ylab("Propensity to Vote") + xlab("") +
    theme(text = element_text(size=18,family = "serif")) + xlab("Coalition Government Lottery") + 
    geom_text(aes(label = paste(ifelse(is.na(est),"",round(est, 2)),
                                ifelse(is.na(stars),"",stars)),y=conf.high),nudge_y = .6)  
  ggsave("FigureSM6.pdf", width = 8, height=6)
  

 
  
  # Regression Tables
  m <- list(m1A_het,m1AFE_he, m1ASC_he,m1B_het,m1BFE_he, m1BSC_he)
  N <- unlist(lapply(m, nobs))
  
  # Table 1
  output <- texreg(m, 
         booktabs = T,
         custom.coef.map=
           list("Intercept" = "Baseline Certain Coal. Lottery",
                "case2" = "Uncertain Coal. Lottery",
                "case3" = "Very Uncertain Coal. Lottery"),
         custom.header = list("CDU/CSU" = 1:3, "Greens" = 4:6),
         custom.model.names = c("Model 1", "Model 2","Model 3", 
                                "Model 1", "Model 2","Model 3"),
        custom.gof.rows = list("Respondent Fixed Effects" = c("","X","","","X",""),
                               "Scenario Fixed Effects" = c("","","X","","","X"),
                               "N" = N)
         )
  cat(output, file = "TableSM21.tex")
  

  # Table 2
  output <- texreg(list(m2A_he,m2B_he), 
         booktabs = T,
         custom.coef.map=
           list("Intercept" = "Baseline Certain Coal. Lottery",
                "case2" = "Uncertain Coal. Lottery",
                "case3" = "Very Uncertain Coal. Lottery",
                "risk_pref(3.33,6.67]" = "Middle Risk Preferences",
                "risk_pref(6.67,10]" = "High Risk Preferences",
                "case2:risk_pref(3.33,6.67]"  = "Uncertain X Middle Risk Pref.", 
                "case3:risk_pref(3.33,6.67]"  = "V. Uncertain X Middle Risk Pref.",  
                "case2:risk_pref(6.67,10]"  = "Uncertain X High Risk Pref.", 
                "case3:risk_pref(6.67,10]"  = "V. Uncertain X High Risk Pref." 
                ),
         custom.header = list("CDU/CSU" = 1, "Greens" = 2),
         custom.model.names = c("Model 1", "Model 1"), 
         custom.gof.rows = list("N" = unlist(lapply(list(m2A_he,m2B_he),nobs))))
  cat(output, file = "TableSM22.tex")
  
  
  
  # Table 3
  output <- texreg(list(m3A_he,m3B_he), 
         booktabs = T,
         custom.coef.map=
           list("Intercept" = "Baseline Certain Coal. Lottery",
                "case2" = "Uncertain Coal. Lottery",
                "case3" = "Very Uncertain Coal. Lottery",
                "ptvA(3.33,6.67]" = "Middle PTV Party",
                "ptvA(6.67,10]" = "High PTV Party",
                "ptvB(3.33,6.67]" = "Middle PTV Party",
                "ptvB(6.67,10]" = "High PTV Party",
                "case2:ptvA(3.33,6.67"  = "Uncertain X Middle PTV Party", 
                "case3:ptvA(3.33,6.67"  = "V. Uncertain X Middle PTV Party",  
                "case2:ptvA(6.67,10]"  = "Uncertain X High PTV Party", 
                "case3:ptvA(6.67,10]"  = "V. Uncertain X High PTV Party",
                "case2:ptvB(3.33,6.67"  = "Uncertain X Middle PTV Party", 
                "case3:ptvB(3.33,6.67"  = "V. Uncertain X Middle PTV Party",  
                "case2:ptvB(6.67,10]"  = "Uncertain X High PTV Party", 
                "case3:ptvB(6.67,10]"  = "V. Uncertain X High PTV Party"
           ),
         custom.header = list("CDU/CSU" = 1, "Greens" = 2),
         custom.model.names = c("Model 1", "Model 1"), 
         custom.gof.rows = list("N" = unlist(lapply(list(m3A_he,m3B_he),nobs))))
  cat(output, file = "TableSM23.tex")
  
  
# Robustness: Interacting with pre-treatment expectations
  
  # Expectation difference averaged over all scenarios
  
  summary(m2A <- lm((vignetteA) ~ case*exp_diff_union_cut, as.data.frame(dta_long)))
  summary(m2B <- lm((vignetteB) ~ case*exp_diff_greens_cut, as.data.frame(dta_long)))
  
  
  output <- texreg(list(m2A,m2B), 
         booktabs = T,
         custom.coef.map=
           list("Intercept" = "Baseline Certain Coal. Lottery",
                "case2" = "Uncertain Coal. Lottery",
                "case3" = "Very Uncertain Coal. Lottery",
                "exp_diff_union_cut(0.757,1.05]" = "Middle Exp. Diff.",
                "exp_diff_union_cut(1.05,1.34]" = "High Exp. Diff.",
                
                "exp_diff_greens_cut(0.759,1.05]" = "Middle Exp. Diff.",
                "exp_diff_greens_cut(1.05,1.35]" = "High Exp. Diff.",
                
                "case2:exp_diff_union_cut(0.757,1.05]"  = "Uncertain X Middle Exp. Diff.", 
                "case3:exp_diff_union_cut(0.757,1.05]"  = "V. Uncertain X Middle Exp. Diff.", 
                "case2:exp_diff_union_cut(1.05,1.34]"  = "Uncertain X High Exp. Diff.", 
                "case3:exp_diff_union_cut(1.05,1.34]"  = "V. Uncertain X High Exp. Diff.", 
                
                "case2:exp_diff_greens_cut(0.759,1.05]"  = "Uncertain X Middle Exp. Diff.", 
                "case3:exp_diff_greens_cut(0.759,1.05]"  = "V. Uncertain X Middle Exp. Diff.", 
                "case2:exp_diff_greens_cut(1.05,1.35]"  = "Uncertain X High Exp. Diff.", 
                "case3:exp_diff_greens_cut(1.05,1.35]"  = "V. Uncertain X High Exp. Diff." 
           ),
         custom.header = list("CDU/CSU" = 1, "Greens" = 2),
         custom.model.names = c("Model 1", "Model 1"), 
         custom.gof.rows = list("N" = unlist(lapply(list(m2A,m2B),nobs))))
  cat(output, file = "TableSM24.tex")
  
  m2A_vcov <- vcovCL(m2A, type="HC1", cluster=dta_long$id)
  m2B_vcov <- vcovCL(m2B, type="HC1", cluster=dta_long$id)
  
  m2A_he <- coeftest(m2A, m2A_vcov)
  m2B_he <-  coeftest(m2B, m2B_vcov)
  
  dfA <- ggpredict(m2A, terms = c("case","exp_diff_union_cut")) %>% as.data.frame() %>% mutate(party = "CDU/CSU")
  dfB <- ggpredict(m2B, terms = c("case","exp_diff_greens_cut")) %>% as.data.frame() %>% mutate(party = "Greens")
  df <- bind_rows(dfA,dfB)
  
  
  
  levels(df$x) <- c("Certain", "Uncertain","V. Uncertain")
  
  df$group <- factor(c("Low Exp. Diff", "Middle Exp. Diff",
                       "High Exp. Diff"), levels=c("Low Exp. Diff", "Middle Exp. Diff",
                                                   "High Exp. Diff"))
  
  # Calcualte Treatment Effects
  ateA_case2_expdifflow <- coef(m2A)["case2"]
  ateA_case2_expdiffmid <- coef(m2A)["case2"] + coef(m2A)["case2:exp_diff_union_cut(0.757,1.05]"]
  ateA_case2_expdiffhgh <- coef(m2A)["case2"] + coef(m2A)["case2:exp_diff_union_cut(1.05,1.34]"]
  
  ateA_case3_expdifflow <- coef(m2A)["case3"]
  ateA_case3_expdiffmid <- coef(m2A)["case3"] + coef(m2A)["case3:exp_diff_union_cut(0.757,1.05]"]
  ateA_case3_expdiffhgh <- coef(m2A)["case3"] + coef(m2A)["case3:exp_diff_union_cut(1.05,1.34]"]
  
  seA_case2_expdifflow <- sqrt(m2A_vcov["case2","case2"])
  seA_case2_expdiffmid <- sqrt(m2A_vcov["case2","case2"] + 
                                 m2A_vcov["case2:exp_diff_union_cut(0.757,1.05]","case2:exp_diff_union_cut(0.757,1.05]"] +
                                 2 * m2A_vcov["case2","case2:exp_diff_union_cut(0.757,1.05]"])
  
  seA_case2_expdiffhgh <- sqrt(m2A_vcov["case2","case2"] + 
                                 m2A_vcov["case2:exp_diff_union_cut(1.05,1.34]","case2:exp_diff_union_cut(1.05,1.34]"] +
                                 2 * m2A_vcov["case2","case2:exp_diff_union_cut(1.05,1.34]"])
  
  seA_case3_expdifflow <- sqrt(m2A_vcov["case3","case3"])
  seA_case3_expdiffmid <- sqrt(m2A_vcov["case3","case3"] + 
                                 m2A_vcov["case3:exp_diff_union_cut(0.757,1.05]","case3:exp_diff_union_cut(0.757,1.05]"] +
                                 2 * m2A_vcov["case3","case3:exp_diff_union_cut(0.757,1.05]"])
  seA_case3_expdiffhgh <- sqrt(m2A_vcov["case3","case3"] + 
                                 m2A_vcov["case3:exp_diff_union_cut(1.05,1.34]","case3:exp_diff_union_cut(1.05,1.34]"] +
                                 2 * m2A_vcov["case3","case3:exp_diff_union_cut(1.05,1.34]"])
  
  # Calcualte Treatment Effects
  ateB_case2_expdifflow <- coef(m2B)["case2"]
  ateB_case2_expdiffmid <- coef(m2B)["case2"] + coef(m2B)["case2:exp_diff_greens_cut(0.759,1.05]"]
  ateB_case2_expdiffhgh <- coef(m2B)["case2"] + coef(m2B)["case2:exp_diff_greens_cut(1.05,1.35]"]
  
  ateB_case3_expdifflow <- coef(m2B)["case3"]
  ateB_case3_expdiffmid <- coef(m2B)["case3"] + coef(m2B)["case3:exp_diff_greens_cut(0.759,1.05]"]
  ateB_case3_expdiffhgh <- coef(m2B)["case3"] + coef(m2B)["case3:exp_diff_greens_cut(1.05,1.35]"]
  
  seB_case2_expdifflow <- sqrt(m2B_vcov["case2","case2"])
  seB_case2_expdiffmid <- sqrt(m2B_vcov["case2","case2"] + 
                                 m2B_vcov["case2:exp_diff_greens_cut(0.759,1.05]","case2:exp_diff_greens_cut(0.759,1.05]"] +
                                 2 * m2B_vcov["case2","case2:exp_diff_greens_cut(0.759,1.05]"])
  
  seB_case2_expdiffhgh <- sqrt(m2B_vcov["case2","case2"] + 
                                 m2B_vcov["case2:exp_diff_greens_cut(1.05,1.35]","case2:exp_diff_greens_cut(1.05,1.35]"] +
                                 2 * m2B_vcov["case2","case2:exp_diff_greens_cut(1.05,1.35]"])
  
  seB_case3_expdifflow <- sqrt(m2B_vcov["case3","case3"])
  seB_case3_expdiffmid <- sqrt(m2B_vcov["case3","case3"] + 
                                 m2B_vcov["case3:exp_diff_greens_cut(0.759,1.05]","case3:exp_diff_greens_cut(0.759,1.05]"] +
                                 2 * m2B_vcov["case3","case3:exp_diff_greens_cut(0.759,1.05]"])
  seB_case3_expdiffhgh <- sqrt(m2B_vcov["case3","case3"] + 
                                 m2B_vcov["case3:exp_diff_greens_cut(1.05,1.35]","case3:exp_diff_greens_cut(1.05,1.35]"] +
                                 2 * m2B_vcov["case3","case3:exp_diff_greens_cut(1.05,1.35]"])
  df_stars <- data.frame(
    est = c(ateA_case2_expdifflow,ateA_case2_expdiffmid,ateA_case2_expdiffhgh,
            ateA_case3_expdifflow,ateA_case3_expdiffmid,ateA_case3_expdiffhgh,
            ateB_case2_expdifflow,ateB_case2_expdiffmid,ateB_case2_expdiffhgh,
            ateB_case3_expdifflow,ateB_case3_expdiffmid,ateB_case3_expdiffhgh),
    se = c(seA_case2_expdifflow,seA_case2_expdiffmid,seA_case2_expdiffhgh,
           seA_case3_expdifflow,seA_case3_expdiffmid,seA_case3_expdiffhgh,
           seB_case2_expdifflow,seB_case2_expdiffmid,seB_case2_expdiffhgh,
           seB_case3_expdifflow,seB_case3_expdiffmid,seB_case3_expdiffhgh),
    group = factor(c("Low Exp. Diff", "Middle Exp. Diff",
                     "High Exp. Diff"), levels=c("Low Exp. Diff", "Middle Exp. Diff",
                                                 "High Exp. Diff")),
    party = c(rep("CDU/CSU",6),rep("Greens",6)),
    x=c(rep("Uncertain",3),rep("V. Uncertain",3))
  )  %>%
    mutate(pval = pnorm(abs(est/se),lower.tail = F)*2 )%>%
    mutate(stars = case_when(
      pval < 0.001 ~ "***",
      pval < 0.01 ~ "**",
      pval < 0.05 ~ "*",
      TRUE ~ ""
    )) 
  
  df <- left_join(df, df_stars, by=c("group","party","x" )) %>%
    as.data.frame()
  
  ggplot(df, aes(x=x,y=predicted,ymin=conf.low,ymax=conf.high)) + 
    geom_errorbar(width=.05) + theme_minimal() +
    geom_point() + 
    facet_grid(group ~ party) +
    geom_bar(stat = "identity", alpha=0.3) +
    ylab("Propensity to Vote") + xlab("Coalition Lottery") +
    theme(text = element_text(size=18)) + 
    geom_text(aes(label = paste(ifelse(is.na(est),"",round(est, 2)),
                                ifelse(is.na(stars),"",stars)),y=conf.high),nudge_y = .6)
  ggsave("FigureSM7.pdf", width = 8, height=6)

  
  # Robustness: Anti-AfD Effect? Excluding AfD vigniettes
  
  dta_long_without_afd_A <- dplyr::filter(dta_long, !grepl('AfD', X_Coal_Text))

  # Basic model
  summary(m1A <- lm(vignetteA ~ case, data = dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]))
  m1A_vcov <- vcovCL(m1A, type="HC1", cluster=dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]$id)
  m1AB_he <- coeftest(m1A,m1A_vcov)
  
  # Fixed Effects
  summary(m1AFE <- lm(vignetteA ~ case + as.factor(id), data = dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]))
  m1AFE_vcov <- vcovCL(m1AFE, type="HC1", cluster=dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]$id)
  m1AFE_he <- coeftest(m1AFE,m1AFE_vcov)
  
  # Scenarios
  summary(m1ASC <- lm(vignetteA ~ case + X_Coal_Text + Y_Coal_Text + Z_Coal_Text, data = dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]))
  m1ASC_vcov <- vcovCL(m1ASC, type="HC1", cluster=dta_long_without_afd_A[which(dta_long_without_afd_A$exp_diff_union_cut == "(0.463,0.757]"),]$id)
  m1ASC_he <- coeftest(m1ASC,m1ASC_vcov)
  
  # Table 4
  output <- texreg(list(m1AB_he,m1AFE_he, m1ASC_he), 
         booktabs = T,
         custom.coef.map=
           list("Intercept" = "Baseline Certain Coal. Lottery",
                "case2" = "Uncertain Coal. Lottery",
                "case3" = "Very Uncertain Coal. Lottery"),
         custom.header = list("CDU/CSU" = 1:3),
         custom.model.names = c("Model 1", "Model 2","Model 3"),
         custom.gof.rows = list("Respondent Fixed Effects" = c("","X",""),
                                "Scenario Fixed Effects" = c("","","X"),
                                "N" = unlist(lapply(list(m1AB_he,m1AFE_he, m1ASC_he),nobs))))
  cat(output, file = "TableSM25.tex")
  